
> # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> # Collective victimhood beliefs and conflict-related attitudes: A meta-analysis
> # Marko Kljajic, Nadav Shelef, and Ethan vanderWilden
> 
> # Last Replicated: August 6, 2025
> 
> # Analysis file
> # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
> 
> ############ 1. Load in packages, data, helper functions ############
> # List of packages
> pkg = c("tidyverse", "patchwork", "metafor", "kableExtra", "rmeta", "reshape2", "metaSEM")

> # Checks if they are installed, install if not
> if (length(setdiff(pkg, rownames(installed.packages()))) > 0) {install.packages(setdiff(pkg, rownames(installed.packages()))) }

> # Load + clear environment
> suppressMessages(suppressWarnings(lapply(pkg, library, character.only = TRUE)))
[[1]]
 [1] "metaSEM"    "OpenMx"     "reshape2"   "rmeta"      "kableExtra" "metafor"    "numDeriv"   "metadat"    "Matrix"     "compute.es" "patchwork"  "lubridate"  "forcats"   
[14] "stringr"    "dplyr"      "purrr"      "readr"      "tidyr"      "tibble"     "ggplot2"    "tidyverse"  "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
[27] "methods"    "base"      

[[2]]
 [1] "metaSEM"    "OpenMx"     "reshape2"   "rmeta"      "kableExtra" "metafor"    "numDeriv"   "metadat"    "Matrix"     "compute.es" "patchwork"  "lubridate"  "forcats"   
[14] "stringr"    "dplyr"      "purrr"      "readr"      "tidyr"      "tibble"     "ggplot2"    "tidyverse"  "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
[27] "methods"    "base"      

[[3]]
 [1] "metaSEM"    "OpenMx"     "reshape2"   "rmeta"      "kableExtra" "metafor"    "numDeriv"   "metadat"    "Matrix"     "compute.es" "patchwork"  "lubridate"  "forcats"   
[14] "stringr"    "dplyr"      "purrr"      "readr"      "tidyr"      "tibble"     "ggplot2"    "tidyverse"  "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
[27] "methods"    "base"      

[[4]]
 [1] "metaSEM"    "OpenMx"     "reshape2"   "rmeta"      "kableExtra" "metafor"    "numDeriv"   "metadat"    "Matrix"     "compute.es" "patchwork"  "lubridate"  "forcats"   
[14] "stringr"    "dplyr"      "purrr"      "readr"      "tidyr"      "tibble"     "ggplot2"    "tidyverse"  "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
[27] "methods"    "base"      

[[5]]
 [1] "metaSEM"    "OpenMx"     "reshape2"   "rmeta"      "kableExtra" "metafor"    "numDeriv"   "metadat"    "Matrix"     "compute.es" "patchwork"  "lubridate"  "forcats"   
[14] "stringr"    "dplyr"      "purrr"      "readr"      "tidyr"      "tibble"     "ggplot2"    "tidyverse"  "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
[27] "methods"    "base"      

[[6]]
 [1] "metaSEM"    "OpenMx"     "reshape2"   "rmeta"      "kableExtra" "metafor"    "numDeriv"   "metadat"    "Matrix"     "compute.es" "patchwork"  "lubridate"  "forcats"   
[14] "stringr"    "dplyr"      "purrr"      "readr"      "tidyr"      "tibble"     "ggplot2"    "tidyverse"  "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
[27] "methods"    "base"      

[[7]]
 [1] "metaSEM"    "OpenMx"     "reshape2"   "rmeta"      "kableExtra" "metafor"    "numDeriv"   "metadat"    "Matrix"     "compute.es" "patchwork"  "lubridate"  "forcats"   
[14] "stringr"    "dplyr"      "purrr"      "readr"      "tidyr"      "tibble"     "ggplot2"    "tidyverse"  "stats"      "graphics"   "grDevices"  "utils"      "datasets"  
[27] "methods"    "base"      


> rm(list=ls())

> # set working directory (NOTE: change for local file path)
> setwd("C:/Users/19785/Desktop/CPS_replication")

> # load in the dataset
> data <- read.csv("KSV_dataset.csv")

> #Load in functions from file
> source("helperFunctions.R")

> ### Main analysis excludes 'bundled' treatments ; emotions/cognitive perspective outcomes
> main <- subset(data, data$include.main == 1) # (keep emotions and cognitive frames)

> main <- subset(main, main$outcome %in% c("Hawkishness", "Ingroup Attachment", "Outgroup Exclusion", "Reconciliation"))

> ############ 2. Main Results ####
> 
> ### Table 1: Number of estimates
> kbl(get_n(main, exclude_mechs = T)$ntab, format = 'latex', booktabs = T)

\begin{tabular}[t]{lrrrrr}
\toprule
  & Hawkishness & Reconciliation & Outgroup.Exclusion & Ingroup.Attachment & iv.total\\
\midrule
Generic & 93 & 100 & 24 & 84 & 301\\
Competitive & 28 & 75 & 51 & 34 & 188\\
Exclusive & 17 & 15 & 26 & 21 & 79\\
Inclusive & 31 & 60 & 55 & 31 & 177\\
Total & 169 & 250 & 156 & 170 & 745\\
\bottomrule
\end{tabular}

> ### Figure 3: Full individual estimate plot
> meta_plot(main, sepColor = T) # 12 x 6 pdf

> ### Figure 4 (A-D): Individual Estimate Plots
> meta_plot(subset(main, main$treatment == "Collective victimhood"), includeEst = T) # 6 x 4

> meta_plot(subset(main, main$treatment == "Competitive victimhood"), includeEst = T) # 6 x 4

> meta_plot(subset(main, main$treatment == "Exclusive victimhood"), includeEst = T) # 6 x 4

> meta_plot(subset(main, main$treatment == "Inclusive victimhood"), includeEst = T) # 6 x 4

> ### In-text estimates (written/mentioned)
> # effect sizes: (see column 1)
> results <- fill_metatable3(main, exclude_mechs = T) 

> results
                  All outcomes     Hawkishness  Reconciliation Outgroup Exclusion Ingroup Attachment
All victimhood  0.274 (0.044)*  0.188 (0.053)*  0.254 (0.070)*     0.190 (0.079)*     0.388 (0.048)*
General         0.333 (0.046)*  0.296 (0.049)*  0.320 (0.059)*      0.404 (0.207)     0.380 (0.065)*
Competitive     0.495 (0.082)*   0.209 (0.171)  0.534 (0.118)*     0.444 (0.138)*     0.553 (0.091)*
Exclusive       0.389 (0.096)*   0.237 (0.182)  0.552 (0.112)*      0.234 (0.132)     0.596 (0.203)*
Inclusive      -0.260 (0.070)* -0.301 (0.108)* -0.495 (0.106)*     -0.139 (0.094)      0.108 (0.139)

> # Overall: 0.274(0.044), Generic: 0.333 (0.046), Competitive: 0.495, 0.82,
> #   Exclusive: 0.389 (0.096), Inclusive: -0.26 (0.070)
> 
> # estimates range from -2.49 to 2.49
> round(range(main$d.directional) , 2)
[1] -2.49  2.49

> #combined Non-Inclusive victimhood (estimate of 0.391, SE = 0.040)
> summary(meta3L(y = subset(main, main$treatment != "Inclusive victimhood")$d.directional, 
+                v = subset(main, main$treatment != "Inclusive victimhood")$var.d, 
+                cluster = subset(main, main$treatment != "Inclusive victimhood")$article.id))

Call:
meta3L(y = subset(main, main$treatment != "Inclusive victimhood")$d.directional, 
    v = subset(main, main$treatment != "Inclusive victimhood")$var.d, 
    cluster = subset(main, main$treatment != "Inclusive victimhood")$article.id)

95% confidence intervals: z statistic approximation (robust=FALSE)
Coefficients:
          Estimate Std.Error   lbound   ubound z value  Pr(>|z|)    
Intercept 0.391448  0.039939 0.313170 0.469726  9.8012 < 2.2e-16 ***
Tau2_2    0.184246  0.013846 0.157109 0.211384 13.3069 < 2.2e-16 ***
Tau2_3    0.098223  0.021565 0.055957 0.140490  4.5548 5.244e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Q statistic on the homogeneity of effect sizes: 12009.84
Degrees of freedom of the Q statistic: 567
P value of the Q statistic: 0

Heterogeneity indices (based on the estimated Tau2):
                              Estimate
I2_2 (Typical v: Q statistic)   0.6356
I2_3 (Typical v: Q statistic)   0.3388

Number of studies (or clusters): 95
Number of observed statistics: 568
Number of estimated parameters: 3
Degrees of freedom: 565
-2 log likelihood: 858.6112 
OpenMx status1: 0 ("0" or "1": The optimization is considered fine.
Other values may indicate problems.)

> ############ 3. Combinations of CVB and Conflict attitudes ####
> 
> ### Figure 5: Main coefficient plot
> full_plots(results, df = main, exclude_mechs = T, Bonf = T, range = 1.5) # 10 x 10 pdf

> #Included in the text: What is the p-value of non-significant results?
> round(pnorm(-abs(0.209 /(0.171))), 3) #Competitive-Hawkishness [p = 0.111]
[1] 0.111

> round(pnorm(-abs(0.237 /(0.182))), 3) #Exclusive-Hawkishness [p = 0.096]
[1] 0.096

> round(pnorm(-abs(0.404 /(0.207))), 3) #Generic-Outgroup Exclusion [p = 0.025] ### 025!!
[1] 0.025

> round(pnorm(-abs(0.234 /(0.132))), 3) #Exclusive-Outgroup Exclusion [p = 0.038]
[1] 0.038

> round(pnorm(-abs(-0.139 /(0.094))), 3) #Inclusive-Outgroup Exclusion [p = 0.070]
[1] 0.07

> round(pnorm(-abs(0.108 /(0.139))), 3) #Inclusive-Ingroup Attachment [p = 0.219]
[1] 0.219

> ############ 4. Consistency across Contexts ####
> 
> ### Figure 6: estimates by setting
> subset_plots(subset_estimates(main, "Setting"), separate = T, range = 2.5, ind_plots = T) # 8 x 7 render for each (2 plots)
$uni

$non_uni


> #Included in the text: Proportion of estimates from US or Israel? (37%)
> round(nrow(subset(main, main$Setting == "Israel" | main$Setting == "United States"))/745, 2)
[1] 0.37

> ### Figure 7: estimates by victimhood dimension
> subset_plots(subset_estimates(main, "victimhood.dimension"), range = 1.5, jump = 0.5) #9 x 6

> ### Figure 8: estimates by conflict type (NOTE: reset factor orders)
> conflict <- subset_estimates(main, "conflict.type")

> conflict$setting <- factor(conflict$setting,  levels = c("Genocide/\nMass Violence",  "Colonialism/\nSlavery", "Armed\nConflict", "Discrimination"))

> subset_plots(conflict, range = 1.5, jump = 0.5)  #9 x 6

> ############ 5. Consistency across Methods ####
> ### Figure 9: Experimental and Observational studies
> 
> exp <- fill_metatable3(subset(main, main$Design.binary == "Experimental"), exclude_mechs = T) 

> exp[1, 1] #REPORTED in paper: overall, d = 0.062, se = 0.03
[1] "0.062 (0.030)*"

> obs <- fill_metatable3(subset(main, main$Design.binary == "Observational"), exclude_mechs = T) 

> obs[1, 1] #REPORTED in paper: overall, d = 0.353, se = 0.056
[1] "0.353 (0.056)*"

> comparison_plot(exp, obs, subset(main, Design.binary == "Experimental"),  subset(main, Design.binary == "Observational"),
+              exclude_mechs = T) #12 x 6

> ############ 6. Publication Bias ####
> 
> ### Figure 10: Regression test for publication bias
> Edata <- subset(main, main$Design.binary == "Experimental")

> Odata <- subset(main, main$Design.binary != "Experimental")

> summary(lm(d ~ sd.d, data = main)) # d = 0.367 + 0.846*SE

Call:
lm(formula = d ~ sd.d, data = main)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6183 -0.3291 -0.1137  0.2163  1.9129 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.36692    0.03364   10.91  < 2e-16 ***
sd.d         0.84590    0.19359    4.37 1.42e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4284 on 743 degrees of freedom
Multiple R-squared:  0.02505,	Adjusted R-squared:  0.02374 
F-statistic: 19.09 on 1 and 743 DF,  p-value: 1.422e-05


> summary(lm(d ~ sd.d, data = Edata)) # d = -0.026 + 1.3576*SE

Call:
lm(formula = d ~ sd.d, data = Edata)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.39185 -0.08237 -0.02223  0.06358  0.73252 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.0256     0.0242  -1.058    0.291    
sd.d          1.3576     0.1182  11.487   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1605 on 187 degrees of freedom
Multiple R-squared:  0.4137,	Adjusted R-squared:  0.4106 
F-statistic: 131.9 on 1 and 187 DF,  p-value: < 2.2e-16


> summary(lm(d ~ sd.d, data = Odata)) # d = 0.418 + 1.196*SE

Call:
lm(formula = d ~ sd.d, data = Odata)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.77250 -0.34630 -0.08794  0.28219  1.77327 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.41845    0.04204   9.954  < 2e-16 ***
sd.d         1.19571    0.25969   4.604 5.14e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4411 on 554 degrees of freedom
Multiple R-squared:  0.03686,	Adjusted R-squared:  0.03512 
F-statistic:  21.2 on 1 and 554 DF,  p-value: 5.135e-06


> ggplot() +
+   geom_point(aes(x = Edata$sd.d, y = Edata$d), alpha = 0.4, color = "darkblue") +
+   geom_point(aes(x = Odata$sd.d, y = Odata$d), alpha = 0.4, color = "darkred") +
+   
+   geom_smooth(aes(x = data$sd.d, y = data$d), method = "lm", color = "black", lwd = 1.5, se = F) +
+   geom_smooth(aes(x = Edata$sd.d, y = Edata$d), method = "lm", color = "darkblue", lwd = 1.5, se = F) +
+   geom_smooth(aes(x = Odata$sd.d, y = Odata$d), method = "lm", color = "red", lwd = 1.5, se = F) +
+   
+   
+   theme_bw() + labs(y = "Effect size magnitude", x = "Estimate SE (converted for effect size)", 
+                     subtitle = "Overall (green): d = 0.367 + 0.846*SE\nExperimental (blue): d = -0.026 + 1.358*SE\n
+                     Observational (red): d = 0.418 + 1.196*SE") +
+   scale_y_continuous(limits = c(0, 2.5)) +
+   theme(text = element_text(size = 15), plot.subtitle = element_blank()) # 12 x 6

> ############ 7. Emotions and Cognitive Perspectives ####
> 
> #get dataset with Emotions and Cognitive Perspectives (but still excluding bundled treatments)
> main_wECP <- subset(data, data$include.main == 1) # (keep emotions and cognitive frames)

> ### Figure 11: estimates on cognitive frames and emotions
> full_plots(fill_metatable3(main_wECP), df = data, only_mech = T, range = 2) # 10 x 4
